/* BoundaryConditions.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
#if !defined(_BoundaryConditions_)
#define _BoundaryConditions_
#include <src/miscellaneous/OS_basic.h>
#include <src/data/ParticleSet.h>
#include <src/data/OS_Multi_Array.h>


struct DirichletParticle : public particle<OS_double>
	,	boost::shiftable<DirichletParticle, std::size_t >     // <<, >>
{
	typedef Dirichlet							Boundary_type;
	particle<OS_double>& operator()(const DirichletParticle& from) {
		OS_DYNAMIC_CHECK(false,"to implement");return *this;
	};

	DirichletParticle(OS_size nb) : particle<OS_double>(nb){};
	DirichletParticle(const particle<OS_double> & from) : particle<OS_double>(from){};
	DirichletParticle(){};
	inline OS_double & operator[](indice_type n) {
		if (n>=0 && n<(indice_type) size())
			return *(_Myfirst + n);
		zero_real = 0; // problem
		return zero_real;
	}
	inline const OS_double & operator[](indice_type n) const {

		if (n>=0 && n<(indice_type)_size)
			return *(_Myfirst + n);
		return const_zero_real;
	};
	inline DirichletParticle& operator <<= (std::size_t how) {
		if (size() <2) return *this;
		if (size() > 1)
			memcpy(_Myfirst, _Myfirst+1,(size()-1)*sizeof(base_value_type));
		(*this)[size()-1]= 0;
		return *this;
	};
	inline DirichletParticle& operator >>= (std::size_t how) {
		if (size() <2) return *this;
		if (size() > 1)
			memcpy(_Myfirst+1, _Myfirst,(size()-1)*sizeof(base_value_type));
		(*this)[0]= 0;
		return *this;
	};
	OS_double zero;
};
template<class T = p_array_2::value_type >
struct DirichletParticle_list : 
	public particles_list<T>
	,	boost::shiftable<DirichletParticle_list<T>, std::size_t >     // <<, >>
{
	typedef Dirichlet							Boundary_type;
	typedef particles_list<T>					Base;
	typedef DirichletParticle_list				type;

	type(const indice_list_type from) : Base(from),zero_particle_(dim()) {
		zero_particle_.assign(0.);
	};
	type(const Base & from) : Base(from),zero_particle_(from.size()) {
	};
	type(OS_size nb, OS_size dim) : Base(nb,dim),zero_particle_(nb) {
	};
	type(const type & from) : Base(from), zero_particle_(from.size())
	{
	};
	type(){};

	inline reference operator[](indice_type n) {
		if (n>=0 && n<size())
			return Base::operator[](n);
			zero_particle_.assign(0.);
			return zero_particle_;
	}
	inline const_reference operator[](indice_type n) const {
		if (n>=0 && n<size())
			return Base::operator[](n);
		zero_particle_.assign(0.);
		return zero_particle_;
	};
	inline type& operator <<= (std::size_t how) {
		OS_DYNAMIC_CHECK(false,"to implement");return *this;
	};
	inline type& operator >>= (std::size_t how) {
		OS_DYNAMIC_CHECK(false,"to implement");return *this;
	};
private:
	mutable value_type zero_particle_;
};

struct PeriodicParticle : public particle<OS_double>
	,boost::shiftable<PeriodicParticle, std::size_t >     // <<, >>
{
	typedef Periodic Boundary_type;
	PeriodicParticle(OS_size nb) : particle<OS_double>(nb){};
	PeriodicParticle(const particle<OS_double> & from) : particle<OS_double>(from){};
	PeriodicParticle(const PeriodicParticle & from) : particle<OS_double>(from){};
	PeriodicParticle(){};
	inline OS_double & operator[](indice_type n) {
			return *(_Myfirst + jmmath::mod<indice_type,size_type>(n,size()));
	}
	inline OS_double & operator[](indice_type n) const {
			return *(_Myfirst + jmmath::mod<indice_type,size_type>(n,size()));
	}
	inline PeriodicParticle& operator <<= (std::size_t how) {
		if (size() <2) return *this;
		value_type local = (*this)[0];
		if (size() > 1)
			memcpy(_Myfirst, _Myfirst+1,(size()-1)*sizeof(base_value_type));
		(*this)[size()-1]= local;
		return *this;
	};
	inline PeriodicParticle& operator >>= (std::size_t how) {
		if (size() <2) return *this;
		value_type local = (*this)[size()-1];
		if (size() > 1)
			memcpy(_Myfirst+1, _Myfirst,(size()-1)*sizeof(base_value_type));
		(*this)[0]= local;
		return *this;
	};
};

template<class T >
struct PeriodicParticle_list : 
	public particles_list<T>
	,boost::shiftable<PeriodicParticle_list<T>, std::size_t >     // <<, >>
{
	typedef Periodic							Boundary_type;
	typedef particles_list<T>					Base;
	typedef PeriodicParticle_list				type;

	type(OS_size nb, OS_size dim) : Base(nb,dim) {
	};
	type(const indice_list_type & from ) : Base(from) {
	};
	type(const Base & from) : Base(from){
	};
	type(){};
	inline reference operator[](indice_type n) {
			return Base::operator[](jmmath::mod<indice_type,size_type>(n,size()) );
	}
	inline const_reference operator[](indice_type n) const {
			return Base::operator[](jmmath::mod<indice_type,size_type>(n,size()));
	};
	inline PeriodicParticle_list& operator <<= (std::size_t how) {
		OS_DYNAMIC_CHECK(false,"to implement");return *this;
	};
	inline PeriodicParticle_list& operator >>= (std::size_t how) {
		OS_DYNAMIC_CHECK(false,"to implement");return *this;
	};
};


/*! Neumann Type Condition
	Let \f$u(y)\f$, with \f$y:=(y_1,..,y_N) \in [0,1]^N\f$.	This class is the null Neumann Condition
\f[
    \partial_{y_0} u(0,y_2,..) = \partial_{y_0} u(1,y_2,...) = 0
\f]
	We may take into account more general boundary condition :
\f[
    \partial_{y_0} u(0,y_2,..) = u_{left}(y_2,..,y_N) ;  \partial_{y_0} u(1,y_2,..) = u_{right}(y_2,..,y_N)
\f]
We will implement this later if we actually need it
*/
template<typename T>
class Neumann_Vector : 
	public OS_array< T>
	,boost::shiftable<Neumann_Vector<T>, std::size_t >     // <<, >>
{
public :
	typedef Neumann_Vector								type;
	typedef OS_array									Base;
	typedef Neumann										Boundary_type;

	explicit Neumann_Vector(size_type nb) : Base (nb) {};
	explicit Neumann_Vector(size_type dim1, size_type dim2) : Base(dim1,dim2) {};
	explicit Neumann_Vector(size_type dim1, size_type dim2, size_type dim3) : Base(dim1,dim2,dim3)  {};
	explicit Neumann_Vector(const indice_list_type& indx) : Base(indx){};
	explicit Neumann_Vector(const smart_ptr_type& data, ptr_type first) : Base(data,first){};

	type(const Base & from) : Base(from) {
	};
	type(Base & from) : Base(from) {};
	type(){};
	type(const type& from)	: Base(from) {};
	template <class T>
	type(const OS_data_array_ref<T,Dim> & from) : Base(from) {
	};
	type& operator =(const type & from) {
		Base::operator=(from);
		return *this;
	};

	virtual reference operator[](indice_type n) {
		if (n>=size())
			return *(begin()+size()-1);
		if (n<0)
			return (*begin());
		return *(begin() + n);
	};
	virtual const_reference operator[](indice_type n) const {
		if (n>=size())
			return *(begin()+size()-1);
		if (n<0)
			return (*begin());
		return *(begin() + n);
	};
	Neumann_Vector& operator <<= (std::size_t how) {
		memcpy(first_, first_+capacity(Dim-1),(capacity()-capacity(Dim-1))*sizeof(base_value_type));
		return *this;
	};
	Neumann_Vector& operator >>= (std::size_t how) {
		memcpy(first_+capacity(Dim-1), first_, (capacity()-capacity(Dim-1))*sizeof(base_value_type));
		return *this;
	};
};

template<typename T>
class Dirichlet_Vector : 
	public OS_array<T>
	,boost::shiftable<Dirichlet_Vector<T>, std::size_t >     // <<, >>
{
public :
	typedef Dirichlet_Vector								type;
	typedef OS_array										Base;
	typedef Dirichlet										Boundary_type;

	explicit Dirichlet_Vector(size_type nb) : Base (nb) {};
	explicit Dirichlet_Vector(size_type dim1, size_type dim2) : Base(dim1,dim2) {};
	explicit Dirichlet_Vector(size_type dim1, size_type dim2, size_type dim3) : Base(dim1,dim2,dim3)  {};
	explicit Dirichlet_Vector(const indice_list_type& indx) : Base(indx){};
	explicit Dirichlet_Vector(const smart_ptr_type& data, ptr_type first) : Base(data,first){};

	type(const Base & from) : Base(from) {};
	type(Base & from) : Base(from) {};
	type(){};
	type(const type& from)	: Base(from) {};
	template <class T>
	type(const OS_data_array_ref<T,Dim> & from) : Base(from) {};
	type& operator =(const type & from) {
		Base::operator=(from);
		return *this;
	};

	template <std::size_t D> reference get_reference () const {
				zero_.resize( (*begin()).size_list()); zero_.assign(0.);	return reference(zero_.data_,zero_.first_);
		};

	template <>	reference get_reference<1> () const {
				zero_ = 0; return zero_;
		};

	virtual reference operator[](indice_type n) {
		if (n <0 || n >= size()) {
			return get_reference<Dim>();
		};
		return *(begin()+n);
	};
	virtual const_reference operator[](indice_type n) const {
		if (n <0 || n >= size()) {
			return get_reference<Dim>();
		};
		return *(begin()+n);
	};
	inline type& operator <<= (std::size_t how) {
		if (size() == 0) return *this;
		if (capacity() > 1)
			memcpy(first_, first_+capacity(Dim-1),(capacity()-capacity(Dim-1))*sizeof(base_value_type));
		if (Dim == 1)
			*(first_+size()-1) = 0.;
		else 
		{
			((*this)[size()-1])= get_reference<Dim>();
		}
		return *this;
	};
	inline type& operator >>= (std::size_t how) {
		if (size() <2) return *this;
		if (capacity() > 1)
			memcpy(first_+capacity(Dim-1), first_,(capacity()-capacity(Dim-1))*sizeof(base_value_type));
		if (Dim == 1)
			*(first_) = 0.;
		else
			((*this)[0]) = get_reference<Dim>();
		return *this;
	};

	mutable value_type zero_;
};

template<typename T>
class Unbounded_Vector : 
	public OS_array<T>
	,boost::shiftable<Unbounded_Vector<T>, std::size_t >     // <<, >>
{
public :
	typedef Unbounded_Vector								type;
	typedef OS_array										Base;
	typedef Unbounded										Boundary_type;

	explicit Unbounded_Vector(size_type nb) : Base (nb) {};
	explicit Unbounded_Vector(size_type dim1, size_type dim2) : Base(dim1,dim2) {};
 	explicit Unbounded_Vector(size_type dim1, size_type dim2, size_type dim3) : Base(dim1,dim2,dim3)  {};
	explicit Unbounded_Vector(const indice_list_type& indx) : Base(indx){};
	explicit Unbounded_Vector(const smart_ptr_type& data, ptr_type first) : Base(data,first){};

	Unbounded_Vector(const OS_array< T> & from) : Base(from) {};
	Unbounded_Vector(){};
	Unbounded_Vector(const Unbounded_Vector& from)	: Base(from) {};

	template <class T>
	Unbounded_Vector(const OS_data_array_ref<T,Dim> & from) : Base(from) {};

	template <std::size_t D> reference get_reference () const {
				infty_ = *begin(); infty_.assign(OS_MAX_REAL);	return reference(infty_.data_,infty_.first_);
		};

	template <>	reference get_reference<1> () const {
				infty_ = OS_MAX_REAL; return infty_;
		};

	virtual reference operator[](indice_type n) {
		if (n <0 || n >= (indice_type) size()) {
			return get_reference<Dim>();
		};
		return *(begin()+n);
	};
	virtual const_reference operator[](indice_type n) const {
		if (n <0 || n >= (indice_type) size()) {
			return get_reference<Dim>();
		};
		return *(begin()+n);
	};
	inline Unbounded_Vector& operator <<= (std::size_t how) {
		if (size() == 0) return *this;
		if (capacity() > 1)
			memcpy(first_, first_+capacity(Dim-1),(capacity()-capacity(Dim-1))*sizeof(base_value_type));
		if (Dim == 1)
			*(first_+size()-1) = OS_MAX_REAL;
		else
			OS_Array_assign<value_type>()((*this)[size()-1],OS_MAX_REAL);
		return *this;
	};
	inline Unbounded_Vector& operator >>= (std::size_t how) {
		if (size() <2) return *this;
		if (capacity() > 1)
			memcpy(first_+capacity(Dim-1), first_,(capacity()-capacity(Dim-1))*sizeof(base_value_type));
		if (Dim == 1)
			*(first_) = OS_MAX_REAL;
		else
			OS_Array_assign<value_type>()((*this)[0],OS_MAX_REAL);
		return *this;
	};

	mutable value_type infty_;
};





template<typename T>
class Mod_Vector : 
	public OS_array< T>
	,boost::shiftable<Mod_Vector<T>, std::size_t >     // <<, >>
{
public :
	typedef Mod_Vector										type;
	typedef OS_array										Base;
	typedef Periodic										Boundary_type;

	explicit Mod_Vector(size_type nb) : Base (nb) {};
	explicit Mod_Vector(size_type dim1, size_type dim2) : Base(dim1,dim2) {};
	explicit Mod_Vector(size_type dim1, size_type dim2, size_type dim3) : Base(dim1,dim2,dim3)  {};
	explicit Mod_Vector(const indice_list_type& indx) : Base(indx){};
	explicit Mod_Vector(const smart_ptr_type& data, ptr_type first) : Base(data,first){};

	Mod_Vector(const Base & from) : Base(from) {};
	Mod_Vector(){};
	Mod_Vector(const Mod_Vector& from)	: Base(from) {};
	template <class T>
	Mod_Vector(const OS_data_array_ref<T,Dim> & from) : Base(from) {};

	virtual typename reference operator[](indice_type n) {
		if (size() == 0) return *begin();
		return *(begin() + jmmath::mod<indice_type,std::size_t>(n,size()));
	};
	virtual typename const_reference operator[](indice_type n) const {
		if (size() == 0) return *begin();
		return *(begin() + jmmath::mod<indice_type,std::size_t>(n,size()));
	};
	Mod_Vector& operator <<= (std::size_t how) {
		if (size() <2) return *this;
		if (Dim == 1) {
			value_type local = *first_;
			memcpy(first_, first_+1,(capacity()-capacity(Dim-1))*sizeof(base_value_type));
			*(first_+size()-1)= local;
		}
		else {
			value_type local = (*this)[0];
			memcpy(first_, first_+capacity(Dim-1),(capacity()-capacity(Dim-1))*sizeof(base_value_type));
			(*this)[size()-1] = local;
		};
		return *this;
	};
	Mod_Vector& operator >>= (std::size_t how) {
		if (size() <2) return *this;
		if (Dim == 1) {
			value_type local = *(first_+size()-1);
			memcpy(first_+1, first_,(size()-1)*sizeof(base_value_type));
			*(first_)= local;
		}
		else {
			value_type local = (*this)[size()-1];
			memcpy(first_+capacity(Dim-1), first_,(capacity()-capacity(Dim-1))*sizeof(base_value_type));
			(*this)[0] = local;
		};
		return *this;
	};
};



template <class T>
struct Dual<T,Periodic> {
	typedef T type;
};
template <class T>
struct Dual<T,Protected> {
	typedef T type;
};
template <class T>
struct Dual<T,Unbounded> {
	typedef typename T::value_type value_type;
	typedef Dirichlet_Vector<value_type> type;
};
template <class T>
struct Dual<T,Dirichlet> {
	typedef typename T::value_type value_type;
	typedef Unbounded_Vector<value_type> type;
};




#endif
